Minimum DFT-weighted norm interpolation of seismic data using FFT

Seismic data interpolation problem can be posed as an inverse problem where from inadequate and incomplete data one attempts to recover the complete band-limited seismic wavefield. The problem is often ill posed due to the factors such as inaccurate knowledge of bandwidth and noise. In this case, regularization can be used to obtain a unique and stable solution. In this abstract, we formulate band-limited data interpolation as a minimum norm least squares type problem. An adaptive DFT-weighted norm regularization term is used to constrain solutions. In particular, the regularization term is updated iteratively using a smoothed version of the peridogram of the estimated data. The technique allows for adaptive incorporation of prior knowledge of the data such as the spectrum support and the shape of the spectrum. Introduction The problem of interpolation/resampling seismic data from its incomplete observations arises in many processing steps that requires a regular sampling. Many methods has been proposed, for examples, prediction error filtering based method (Claerbout, 1992; Splitz, 1991), dip moveout based interpolation (Biondi, 1999), Fourier based interpolation (Cary, 1997; Duijndam et al., 1999; Hindriks et al, 1997). Among those methods, Fourier based reconstructions do not make geological (geophysical) assumpions other than the data to be reconstructed are spatially band-limited. The methods start by posing the interpolation/resampling problem as an inversion problem where from inadequate and incomplete data one attempts to recover Fourier coefficents of complete seismic wavefield. However, the problem is often effectively underdetermined which, as is well known, can be satisfied by many solutions. In this case, a regularized solution can be used where the regularizer (weighting function) serves to impose a particular feature on the solution. The criteria to choose a suitable weighting function have been discussed by several researchers (Cabrera and Thomas, 1991; Duijindam et al., 1999; Hindriks et al., 1997; Sacchi and Ulrych, 1998; Zwartjes and Duijndam, 2000). The essential ideal is to incorperate as much limited prior knowledge as possible. e.g. a regularizer imposed by the Cauchy distribution can be used to obtain an estimation of Fourier transform with sparse distribution of spectral amplitudes in Fourier domain (Sachhi and Ulrych, 1998); a weighting function based on the distance between samples or the area surrounding samples is effective in nonuniform Fourier reconstruction of irregular sampled data along one or two spatial coordinates (Duijindam et al., 1999; Hindriks et al., 1997). In this abstract, a band-limited seimic data interpolation method is developed in space domain. The method is equivalent to a Fourier domain reconstruction when discrete Fourier transform is used to obtain the solution (Liu and Sacchi, 2001). In particular, we have modified the least squares approach to include an adaptive DFT-weighted norm regularization term which incorporates a priori knowledge of energy distribution in wavenumber domain that deals with the ill posed band-limited signal interpolation problem. And, unlike conventional method, the bandwidth of the data is not assumed to be known in our approach. An adaptive frequency weighted norm scheme has been proposed by Cabrera and Parks (1991) to extropolate time series. In their approach, the method of modified periodiom is used to obtain adaptive weights from a previous estimation of the time series. In our approach, the weighting function is updated through iterations using a smoothed version of the periodogram. The smoothing was done by convolving a suitable function which is useful as a mean to reduce to irregularity in the spectrum introduced by missing samples. In addition, we show that the iterative regularization can be done very efficiently using FFT and preconditioned conjugate gradient algorithm. The new method can be applied to seismic data in any domains with one or two spatial coordinates. Finally, examples illustrate effectiveness of the method for 3D real seismic data interpolation. Adaptive weighted regularization We will denote y the lengthM vector of observations and x the lengthN vector of unknowns such that where T is the M N × sampling matrix of the problem. The interpolation problem can, therefore, be posed as an inverse problem where from inadequate and incomplete data y one attemps to recover complete data x . Note that the problem is rank deficient. The uniqueness of solution of the problem can be imposed by defining a regularized solution by solving the problem which is often expressed by where 2 || || ⋅ stands 2 l norm, μ is a specified weighting factor controls the trade off between the data norm and misfit of observations. In this abstract, we have modified the regularization in (2) using DFT-weighted norm, in which case the particular object function to be minimized is y x = T (1) 2 2 2 2 ( ) || || || || wr P J x x x y μ = + − T . (3) 2 2 2 2 2 ( ) || || || || , J x x x y μ = + − L T (2)